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Abstract 

We investigate the free cooling of inelastic rough spheres in the pres- 
ence of Coulomb friction. Depending on the coefficients of normal restitu- 
tion e and Coulomb friction n, we find qualitatively different asymptotic 
states. For nearly complete normal restitution (e close to 1) and large fi, 
friction does not change the cooling properties qualitatively compared to a 
constant coefficient of tangential restitution. In particular, the asymptotic 
state is characterized by a constant ratio of rotational and translational 
energies, both decaying according to Haff's law. However, for small e and 
small /i, the dissipation of rotational energy is suppressed, so that the 
asymptotic state is characterized by constant rotational energy while the 
translational energy continues to decay as predicted by Haff's law. Intro- 
ducing either surface roughness for grazing collisions or cohesion forces 
for collisions with vanishing normal load, causes the rotational energy to 
decay according to Haffs law again in the asymptotic long-time limit with, 
however, an intermediate regime of approximately constant rotational en- 
ergy. 

1 Introduction 

Impact properties of small grains have been measured by several groups ^ 
The experimental data is frequently parametrized using a simple model intro- 
duced by Walton The model involves three parameters; the first one, e 
characterizes the incomplete restitution of the normal component of the rela- 
tive velocity of the contact point, denoted by g. The second one, Coulomb's 
coefficient of friction /i describes the reduction of the tangential component of 
g due to sliding, while the third parameter, (3q accounts for the incomplete 
restitution of the tangential component of g for sticking contacts. All three 
parameters have been measured experimentally for various materials making a 
well calibrated model available for theoretical investigations. 
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2 Binary collisions 



We have previously investigated the free coohng of rough spheres ^ and 
needles Q using the formalism of a Pseudo-Liouville operator. In this paper we 
extend the analysis to include Coulomb friction. We briefly recall the collision 
rules for two spheres of equal diameter a, mass m and moment of inertia /. The 
unit- vector from the center of sphere two (r2) to the center of sphere one (ri) 
is denoted by h := (ri — r2)/|ri — r2|. Center-of-mass velocities and angular 
velocities before collision are denoted by Vi, V2, uJi and u)2- Post-collisional 
quantities are primed. The relative velocity of the contact point is given by 
g = Vi — V2 + X {u)i + UJ2)- The relative velocity after collision is given by 

n- g ^- e{g,n){h- g) with e(g, n) e [0, 1], (1) 

nx g' ^- (3{g,n){nx g) with /3(g, n) e [-1, 1] (2) 

where e{g,n) and f3{g,h) are the coefficients of restitution, which in general 
depend on g and n. We assume e to be constant and allow f3 to depend on 
the angle 7 between g and h in order to account for the different energy loss 
mechanisms of sliding and sticking contacts. The impact angle satisfies 7 € 
tt], so that C0S7 = n • g/\g\ < 0. 

The two constitutive equations plus the conservation laws for linear 

and angular momenta determine the post-collisional velocities 

Vi = Vi + Av, V2 = V2 — Av, 

u)i + Au), u)2 ^ <^2 + Au; (3) 



where 



Av = — (n • Vi2)n - ijn x [vi2 x n + -u;i2 

Alj =^ (n X V12 + X {nx lju)^ (4) 

with V12 = f 1 - V2, a;i2 = a;i + a;2, V = 7/(7) 17(1 + /3(7))/(2(l -|- q)), and 
q :— (4/)/(ma^) (q — 0.4 for homogeneous spheres). 

Sliding contacts are governed by Coulomb friction, giving rise to an impact 
angle dependent coefRcient of tangential restitution 

/3*(7)--l-^(l + e)MC0t7. (5) 

q 

Sliding, however, only occurs for impact angles 7 < 70- For higher values 
one expects sticking which is governed by a constant coefficient of tangential 
restitution —l<Po<l such that (3o is strictly > — 1. In agreement with 
Walton § we assume that either sliding or sticking occurs in any single collision, 
never both. We require (3{j) to be continuous and obtain for the limiting angle 
70 = 7o(e,/3o,M) 

q 1 + Po l 

cot7o = (6) 

' 1 + q 1 + e ^l ^ ^ 
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so that 

/3(7) = min{/3o,-l-^(l + e)^cot7} (7) 

As shown in Fig.(|^), /i ^ corresponds to smooth spheres and fi ^ oo amounts 
to constant tangential restitution. 
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Figure 1: (3{'-f) for different values oi fi; e = 0.5, /3o = 0.5. 



3 Free Cooling of the many particle system 
3.1 Analytical Theory 

We consider a system of N classical particles confined to a 3-dimensional volume 
V interacting through a hard-core potential. The time evolution of a dynamic 
variable A = A{{rk{t),Vk{t),uJk{t)}) is determined by a pseudo-Liouville oper- 
ator for t > 

A{{rk,Vk,u}k}, t) = exp(i/:+<)A({rfe, t>fe, Wfe}, 0). (8) 

The pseudo-Liouville operator C+ consists of two parts £+ = Co + C^. The 
first one, describes the free streaming of particles 

Co^-lY^Vk-Vr,, (9) 
k 
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and the second one, = k^k^i'^+i^^) describes hard-core collisions of two 
particles 

T+{kl) = i{vki ■ f-kiM-Vki ■ rki)S{\rki\ - a){b+ - 1). (10) 

The operator replaces the linear and angular momenta of two particles k 
and I before collision by the corresponding ones after collision, according to eqs. 
(H). Q{x) is the Heaviside step-function, and we have introduced the notation 
fki = f"fc — ri and rki — rki/\rki\- Equation ( p^ ) has a simple interpretation. 
The factor v^i ■ Vki gives the flux of incoming particles. The 8- and 5-functions 
specify the conditions for a collision to take place. A collision between particles 
k and / happens only if the two particles are approaching each other which is 
ensured by Q{—Vki ■ th)- At the instant of a collision the distance between the 
two particles has to vanish expressed by S{\rki \—a)- Finally, (6^; — 1) generates 
the change of linear and angular momenta. 

The ensemble average of a dynamic variable is defined by 



{A)t - / dTp{0)Ait) - / dTpit)A{0) 

(11) 

'[[{drkdvkduJk)p{t)A{0). 

k 

Here p{t) — exp {—iC\_t) p{Q) is the A^-particle distribution function, whose time 
development is governed by the adjoint of the time evolution operator C+. 
Differentiating equation (O) with respect to time we get 



= J dTp{0)^^A{t) = J dTpmC+Ait) 
= J dTp{0)exp{iC+t)iC+A{0) (12) 
= J drp{t)iC+AiO) = {iC+A)t 
We are interested in the average translational and rotational energies per particle 

i 



N ^ 2 



N ^ 2 

i 

as well as the total kinetic energy E ~ En + Ej-ot ■ 

We assume a homogeneous cooling state (HCS) and approximate the TV- 
particle distribution function by a Gaussian 

.«)o<ne(h,|-«)e.p{-(5f|^ + 5;^)} (15) 
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where the product of Hcaviside functions accounts for the excluded vohime. 
The state of the system depends on time only through the average translational 
and rotational kinetic energies. Hence its full time dependence (within the HCS 
approximation) is determined by two coupled differential equations for Ttr(f) 
and Trot (i) 

Q J ^ 

2jt^rot{t) =J^{Erot)t = {lC+E,ot)t (17) 

The expectation values on the right hand side can be calculated for the HCS 
state. We obtain 



2 (1 + if;)(l + ^OS' ^0 + C0S2 7o) 




(l + 2k.cos2 7o)2 

sin7o 



^ arctan(^l + 3j^cot7o^ 
l + 31^cos27o+ 



■ cos 70 



.3/2^0 [ 2r?o(l + |e) 
2g [l + 31^cos2 7o 

7^ / sin7o ^ arctan(^l + |^cot7o) \ 1 
^t-- Vl + it^°^'To yr + ^cos7o // 

where v — IG/'i^Jn /ma'^nog{a) sets the time scale, r]Q = 77 sin 70 = [q{l + 
/3o}/{2(l + g)}]sin7o, g{a) denotes the pair correlation at contact, and no = 
N/V. 

Rotational energy is conserved only in the perfectly smooth case character- 
ized by (/3o = —1) or (^ = 0). Translational energy is conserved only if in 
addition e = 1. The total energy is conserved if both translational and rota- 
tional energies are conserved or in the perfectly rough case [n — 00 /\ (3q = 
with complete normal restitution e = 1. For all other values of the parameters 
e, /3o, and ^, the translational and rotational energies decrease with time. 

3.2 Simulations 

Simulations are performed using an event-driven algorithm where the particles 
follow an undisturbed translational and rotational motion until a collision oc- 
curs. In a collision, the particles' velocities just after contact are computed 



V at 
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using the velocities just before contact as stated in cqs. (||). To accelerate the 
simulations we use the algorithm of Lubachevsky and a linked cell structure, 
which allows us to look for collision partners in the neighborhood of a given 
particle only. 

To obtain a well-defined initial configuration we start the simulation on a 
regular lattice with random velocities chosen from a Boltzmann distribution and 
zero angular velocities. To equilibrate the system we choose e = 1 and fio = —1 
corresponding to perfectly smooth spheres and let the simulation run for 200 
collisions per particle. Then e, /i, and Po a-re switched to their desired values. 

To circumvent the problem of inelastic collapse, i.e. the time between two 
collisions becomes too short to be resolved properly, we use the tc-model [|[: if 
the time between a collision and the preceding one for at least one particle is 
smaller than a critical value tc, the collisions parameters are set to their elastic 
values. 

We are not primarily interested in phenomena like shear and cluster in- 
stabilities, but want to investigate how friction effects the cooling properties 
in the rapid flow regime. Hence we simulate dilute systems and aim at good 
statistics. We perform simulations of 3250 particles with a volume fraction 
p = ^a'^^ — 0.101. For all plots we introduce the dimensionless temperatures 
T = Ttr/Ttr(0) and R = Trot/Ttr(0) and a dimensionless time r = iyy^TUP)t . 
The pair correlation function at contact, g{a) is computed using the Carnahan- 
Starling formula in 3D ||^: 

= • (19) 

All data presented here corresponds to initially non-rotating (-R(O) = 0, T(0) = 
1) homogeneous spheres {q = 0.4). 



3.3 Comparison of analytical theory and numerical simu- 
lations 

The most surprising result for the model with impact angle dependent tangential 
restitution is a transition which separates two phases with different asymptotic 
decays of translational and rotational energies. For large fi and nearly complete 
normal restitution (e close to 1), we observe cooling properties which are very 
similar to those obtained in the model with constant tangential restitution cor- 
responding to the limit of infinite /i. The asymptotic state is characterized by 
a constant ratio of rotational to translational energy, both decaying according 
to Haff's law The time dependence of the translational and rotational en- 
ergies for a typical set of parameters, corresponding to soda lime glass Q, is 
shown in Fig. ^. The numerical solution of eqs. (^8|) is compared to simulations 
and good agreement is found for the whole cooling range. For short times there 
is a linear change of both temperatures, whereas in the asymptotic state both 
temperatures decay like t"^ according to Haff's law. In the asymptotic state 
the ratio of rotational and translational temperature is constant in time. 
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Figure 2: Decay of translational and rotational energy for a set of parameter 
values (e = 0.97, /i = 0.092, /Jq = 0.44) corresponding to Soda lime glass; results 
of the approximate analytical theory are compared to data from simulations. 



For small ^ and small e the rotational energy remains constant in time 
(after an initial increase for small initial i?(0)). The translational energy decays 
according to Half's law. An example is shown in Fig. ||for e = 0.3 and /i = 0.2. 
This, at first surprising result can be explained quite easily: Coulomb's law 
of friction yields only very small friction for small normal loads. So, when the 
spheres lose large amounts of their translational energy (small e) but only a tiny 
bit of their rotational energy (small /i) the system develops towards a state in 
which hardly any more rotational energy is lost because the grains only suffer 
impacts with very small normal load. Hence the asymptotic state resembles that 
of smooth spheres: The system consists of a finite fraction of rotating particles 
at rest. To discuss this state analytically we expand eqs. (|l8|) for small /i, which 
implies 70 close to tt. We set 70 = tt — (5 and expand to leading order in 6 

T^Pi^—^) (20) 
^<5(l+/3o) Tt,r,ot 

8(1 + 9) yiV+iWg 

In leading order, the translational energy is decoupled from the rotational energy 
and equal to the one for smooth spheres. The solution of the second equation 
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Figure 3: Comparison of analytical theory and simulations for a set of parame- 
ters e = 0.3, /i = 0.2, and (3o = 0.5 such that the rotational energy survives and 
only the translational energy decays like (Haff's law). 



is given by 



8{l + l3o)TrS^ 1 



The constant is difficult to evaluate because it depends on the time scale of the 
crossover to the asymptotic regime and on the values of Ttr and Trot on this 
time scale. Analytical theory and simulation agree well in this range of parame- 
ters, too. These two regimes with qualitatively different long time behavior are 
indicated in Fig. ^ as rough and smooth. Why there are two critical values fXi 
and fi2 will be explained in the following discussion. 

To locate the transition between these two phases, we investigate the fol- 
lowing question: For which range of parameters do eqs. allow for a solu- 
tion with a constant ratio of rotational to translational energy? We plug the 
ansatz k = Trot /Ttr = R/T into eqs. and use k — ^^fj^- Introducing 

X := y^l -I- ^ we obtain a function f{x) whose zeros are possible solutions for 
an asymptotic state. 

1 + 2x2 cot^ 70 1 



fix) = 277^ 



1 + x"^ cot^ 70 V 1 + cot^ 70 g2(a;2 — 1) 

^ ^^ 1 - g ( 1 ^ arctan(a;cot 70 ) \ _ 1 - 

q \ 1 + a;^ cot^ 70 x cot 70/2 
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Figure 4: Transition lines ^i(e,/3o) and /^2(e) for jSg = 0.5. While ^2 is inde- 
pendent of Po, /ii has a weak dependence on Po- The intermediate phase lying 
between /ii and /i2 is largest for /3o = 1, and /ii ^ /X2 as /3o ^ 1. For /3o = 1 
the curve for fii lies less than 0.01 below the one for /3o — 0.5. 



In the limit fi 00 (70 — + the equation f{xo) ~ reduces to a quadratic 
one for which exactly one positive zero exists. We obtain k — X + Va^ + I 
where 

2772 4 



X = t;^ —r- + v'—^-V—^ (24) 



in agreement with ||, |Tl|, |T^. For < /i < 00 and /3o 7^ — 1 we get < 
I cot7o| < 00. For X > 1, consider f{x) in the limits x I and x 00. 

fix) ^ -00 (25) 



hm fix) = - 1^ = (1 + ,)2^,2 _ If! (26) 

x^oo cot 7o 2 2 



From eqs. (|25|, |2q ) we see that at least one solution xq for /(.tq) = and hence 
for a constant asymptotic ratio R/T = k = qix'j^ ~ 1) exists if 



M>^^^=:/^.(e) (27) 

This critical value is independent of /3o. To find out if there is more than one 
solution, and if there are solutions even if eq. (|2^) is violated we take a look at 
gix) := xfix). Since a; > 1, xq is a zero of fix) if and only if xq is a zero of gix). 
For a given /3o there are three qualitatively different shapes of gix) depending 
on /i (see Fig ||). For small fx ip < fix), gix) has no zero. For fii < /i < /i2, 
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Figure 5: Plot of the function g{x) whose zeros are possible asymptotic states 
for different values of fi. (e — 0.9, (3o = 0.5). In this case /ii = 0.08507 and 
/X2 = 0.16222. For /i ^ /i2, g{x oo) ±00. 



g{x) has two zeros, and for fi > ^2, g{x) has one zero. When g{x) has two zeros 
only the smaller one serves as an asymptotic ratio, and only if the initial value 
i?(0)/T(0) is smaller than the greater zero. If the initial value is greater than 
the greater zero the system behaves like in the regime where no solution for an 
asymptotic ratio exists, that means the rotational energy survives. 

The critical lines /ii(e) and /i2(e) are shown in Fig. ^ for /3o = 0.5. /i2 is 
given by eq. (|2^) and /ii is evaluated numerically. To conclude, we observe 
three different phases: 1) a rough phase with a constant ratio of translational 
and rotational energies, b) a smooth phase with constant rotational energy 
and c) an intermediate phase where the asymptotic state is determined by the 
initial value of rtr(0)/Tiot(0). The intermediate phase is largest for /3o = 1, and 
Ml as /3o ^ 1. 

For /3o = 0.5 the asymptotic ratio R/T is shown in Fig. ^ as a function 
of e for various /i, 0.01 < fi < 00. For fi > fj,i{e = 0) the asymptotic state is 
characterized by a constant ratio of rotational and translational energies for all 
values of e. When fi is decreased to smaller values, we find an asymptotically 
constant ratio only for sufficiently large e. We have chosen i?(0) = so that 
there is only one transition at fii and we do not observe the intermediate phase, 
in which the asymptotic state depends on the choice of initial condition for R(0). 

In Fig. we show the crossover between the two phases with different 
long time asymptotics for the rotational energy. The coefficients of normal and 
tangential restitution are fixed (e — 0.97, /3o = —0.99) while /i is varied over the 
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Figure 6: Asymptotic ratio R/T as a function of e for different values of /x, 
f3o = 0.5. For large fi there exists a constant ratio for all < e < 1. For small 
fi, however, there is a critical e such that for smaller e the rotational energy 
survives while the translational energy decays like t~'^. 
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Figure 7: Crossover between the two phases with different long time asymptotics 
as a function of for fixed e = 0.97 and Po = -0.99. /xi = 0.08717. 
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full range < /i < oo. For small ^ the translational energy is almost the same 
as that for smooth spheres and the rotational energy is constant for long times. 
For /I > /ii the rotational energy decays, the sooner the larger ^. The decay of 
the translational energy is slowest for /i fi'^ . As /i ^ oo it approaches the 
curve for constant tangential restitution. 

Deviations between the approximate analytical theory and simulations are 
observed in the parameter regime, close to the transition lines. In particular, 
the parameters fj, and e can be chosen such that the analytical theory predicts a 
Haff type decay of the rotational energy, whereas the simulations reveal constant 
R. In Fig. ^ we show results of a simulation for the parameters of cellulose 
acetate spheres as measured by Foerster et al. [0. Looking at single grains in 
the simulation, one finds extremely non-Gaussian states, in the sense that few 
particles rotate with high angular velocities and dominate the average rotational 
energy. In Fig. ^ we plot a histogram of the rotational velocities for a snapshot 
taken aX t — 10^. It reveals clearly that the rotational energy is dominated by 
few particles with high rotational velocities. Snapshots taken at other instants 
of time show that the identity of particles with high rotational velocities is 
conserved. The HCS approximation is not expected to hold in such a state 
which is strongly non-Gaussian and has a nearly log normal distribution of the 
angular velocity. 
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T 

Figure 8: Intermediate regime of parameters (e = 0.87, = 0.25, /3o — 0.43) cor- 
responding to cellulose acetate; theory predicts a constant ratio for translational 
and rotational energies, wheras simulations show a time persistent rotational en- 
ergy. 
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Figure 9: Histogram at t — 10^ of ln(|u;|) for cellulose acetate (e = 0.87, 
H = 0.25, Pq = 0.43) 



4 Modification of Coulomb's law for low load 

The persistence of rotational energy during cooling can be traced back to Coulomb's 
law which predicts vanishing frictional losses for grazing collisions, regardles of 
the magnitude of the relative tangential velocity of the contact point (see Fig. 
|l|). For realistic materials one would expect some residual friction due to sur- 
face roughness. This effect can be modeled crudely by a minimal roughness /3min 
such that /3(7) > /?min- In addition. Coulomb's law, |Ffiic| = • -FioadI, which 
we have used only holds for sufficiently large normal loads. When the normal 
load gets very small, as it happens in the late stages of cooling and in particular 
for small e, cohesion begins to play a role so that there will always be friction 
even for zero-load impacts as discussed by Johnson, Kendall, and Roberts [|l^ . 
They modify Coulomb's law according to: 



fric I 



M l" - -FloadI +Fo + \ 2\h ■ FioadI + (28) 



where Fq = | Trails > and Es denotes the material-specific surface energy. 
Even simpler is the following version discussed as early as 1934 llj] 



iFfricI =M(^l"-^load|+Foj (29) 

with a small positive quantity (Fq) due to cohesion. To estimate the effects of 
cohesion we integrate eq. (E9) over the duration of a collision Ai ^ ac~^/^v~'^^^ 
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as given by the theory of Hertz (c denotes the velocity of sound). We then 
obtain a modified Coulomb law 

m\n X Av\ = fi{m\h ■ Av\ + FoAt) (30) 

which corresponds to an impact angle dependent coefficient of tangential resti- 
tution for sliding contacts 

q I m\g x n\ J 

and generalizes eq. (|^) to include cohesion forces for low impact collisions. The 
most important feature is a finite coefficient of tangential restitution for impact 
angles 7 ~ n/2. 

The question arises, whether the time persistence of the rotational energy 
survives, if finite /3(-|) is taken into account as predicted by surface roughness 
as well as by cohesion forces. In Fig. ^ we show results of simulations with an 
angle dependent coefficient of restitution, as given in eq. (j^), however, with a 
lower bound /3min as predicted by surface roughness as well as cohesion forces. 
The rotational energy shows a plateau for intermediate times and decays asymp- 
totically like t~'^ for long times. The length of the plateau and the onset of the 
decay depend on the value of /?min as expected: The plateau is longer and the 
decay sets in at later times the closer /3min is to —1. 

5 Summary and Outlook 

We have investigated the effects of friction on the cooling properties of granular 
particles. We observe three distinct phases which differ qualitatively in their 
late stage of cooling. In the rough phase cooling is characterized by a constant 
ratio of translational and rotational energies whereas the smooth phase is char- 
acterized by a time persistent rotational energy even for the latest times. These 
two regimes are separated by an intermediate regime in which the late stage of 
cooling can be either smooth or rough depending on the initial conditions. Both 
regimes are also observed in the simulations. In fact, approximate analytical 
theory and simulation agree well within both phases. Close to the intermedi- 
ate regime we find a strongly non-Gaussian angular velocity distribution which 
causes the analytical theory to fail. Deviations between theory and simultaion 
may also be due to finite size effects which are expected to also play a role 
in experiments on granular media, in contrast to experiments on conventional 
systems of statistical mechanics, where finite size effects are usually negligible. 

Any small friction for grazing collisions as generated for example by surface 
roughness causes the rotational energy to decay. However, a plateau survives 
for intermediate time scales, such that the decay of the rotational energy sets in 
at much later times than the decay of the translational temperature. The length 
of the plateau of the rotational energy, and hence the time of decay diverges as 
the roughness goes to zero. 
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Figure 10: Simulations for e = 0.3, /i = 0.2, and /3o = 0.5 for different /3min- 

The values of /3,nin arc shown at the corresponding curves. Any — 1 < /3ini„ < /3o 
causes the rotational energy to decay. As /?min approaches -1 the plateau of the 
rotational energy persists for longer times. /3min = /?o cancels all /U-dependence 
of /? and thus reveals a constant /?. 
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A possible extension of our work are driven systems. A particular way of 
driving - adding random velocity vectors - has been investigated recently with 
simulations and approximate analytical theory by Cafiero et al. [ p^ . They find 
an asymptotic state in which both kinetic energies take on constant values due 
to driving. 
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